//*************************************************************************//
//Programm to compute minimum distance to the x closest firms and the closest
//land border crossing with the US
//(c) Behrens, Bougna, and Brown
//
//First version: 07/05/2014
//This  version: 07/05/2014
//
//Updates:
//No updates for now, first version of the program.
//*************************************************************************//


#include <iostream>
#include <gsl/gsl_sf_bessel.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_poly.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_sf.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_statistics_double.h>
#include <gsl/gsl_sort.h>

#include <vector>


//UTILITY MACROS

#define DEG2RAD(x) (x * M_PI) / 180.0;
#define EARTH_RADIUS 6378.137

using namespace std;

enum LOGIC {FALSE = 0 , TRUE = 1};



//UTILITY FUNCTIONS

//Allocate array

double** alloc_large_array(long rows, long cols)
{
    double **array_to_alloc = new double*[rows];
    
    for (long i = 0; i < rows; i++)
    {
        array_to_alloc[i] = new double[cols];
        
        if (array_to_alloc[i] == NULL)
        {
            fprintf(stderr, "ERROR - Memory could not be allocated, aborting\n");
            return(NULL);
        }
    }
    return(array_to_alloc);
}


double** free_large_array(long rows, double **ptr)
{
    for (long i = 0; i < rows; i++)
    {
        free(ptr[i]);
    }
    
    free(ptr);
    
    return(NULL);
}


//MAIN PROCEDURE


int main (int argc, char *const argv[])
{
    //Main parameters of the program
    
    char ofile_out[1024];	//filename used for output
    
    int t, geosize;
    int number1, number2, number3, number4, number5;
    int bordersize;
    
    clock_t start, finish;
    
    double ilat  = 0, ilon  = 0, jlat  = 0, jlon  = 0, distance = 0, blat, blon;
    
    
    //Read in the command line parameters
    
    start = clock();
    
    if (argc < 9)
    {
        fprintf(stderr, "%s: ERROR - Not all arguments correctly specified.\n", argv[0]);
        return(1);
    }
    
    
    //****************************************
    //STEP [1]: Get the command line arguments
    //****************************************
    
    //The universe of the firms to be used is passed to the program as the first argument argv[1]
    
    //Structure of the file is as follows (9 variables):
    //Line 1    : Number of geographical units
    //Line x > 1: locID latitude longitude primarynaics employment (we won't use the latter two!!)
    
    FILE *ifile_geodata;
    double **geospace, **outspace;
    
    if ((ifile_geodata = fopen(argv[1], "r")) == NULL)
    {
        fprintf(stderr, "%s: ERROR - File %s not found, aborting.\n", argv[0], argv[1]);
        return(1);
    }
    
    //Output is written to the file specified as the second argument argv[2]
    strcpy(ofile_out, argv[2]);
    
    //The land border crossing file is passed as the third argument argv[3]
    
    FILE *ifile_lbc;
    
    if ((ifile_lbc = fopen(argv[3], "r")) == NULL)
    {
        fprintf(stderr, "%s: ERROR - File %s not found, aborting.\n", argv[0], argv[2]);
        return(1);
    }
    
    //The number of nearest firms to use is the next five arguments argv[4]-argv[8]
    if ((number1 = atoi(argv[4])) < 1)
    {
        fprintf(stderr, "%s: ERROR - invalid number of firms for computing min dist, aborting.\n", argv[0]);
        return(2);
    }
    
    if ((number2 = atoi(argv[5])) < 1)
    {
        fprintf(stderr, "%s: ERROR - invalid number of firms for computing min dist, aborting.\n", argv[0]);
        return(2);
    }
    
    if ((number3 = atoi(argv[6])) < 1)
    {
        fprintf(stderr, "%s: ERROR - invalid number of firms for computing min dist, aborting.\n", argv[0]);
        return(2);
    }
    
    if ((number4 = atoi(argv[7])) < 1)
    {
        fprintf(stderr, "%s: ERROR - invalid number of firms for computing min dist, aborting.\n", argv[0]);
        return(2);
    }
    
    if ((number5 = atoi(argv[8])) < 1)
    {
        fprintf(stderr, "%s: ERROR - invalid number of firms for computing min dist, aborting.\n", argv[0]);
        return(2);
    }
    
    
    //**************************************
    //STEP [2]: Read in and prepare the data
    //**************************************
    
    
    //Read in data for the universe of firm locations (*geodata file)
    
    fscanf(ifile_geodata, "%d", &geosize);
    geospace = alloc_large_array(geosize, 5);
    
    //This holds the plant identifier, the five mindistances, and the min distance to the US
    outspace = alloc_large_array(geosize, 7);
    
    fprintf(stdout, "%s: Processing the file with the plants...\n", argv[0]);
    
    fprintf(stdout, "%s: Number of plants to be processed = %d\n", argv[0], geosize);
    fprintf(stdout, "%s: Computing across %d %d %d %d %d nearest plants\n", argv[0], number1, number2, number3,
            number4, number5);
    
    t = 0;
    
    while (fscanf(ifile_geodata, "%lf %lf %lf %lf %lf", &geospace[t][0], &geospace[t][1], &geospace[t][2], &geospace[t][3], &geospace[t][4]) != EOF)
    {
        //Convert the latitude and the longitude into radians for distance computations
        geospace[t][1] = DEG2RAD(geospace[t][1]);
        geospace[t][2] = DEG2RAD(geospace[t][2]);
        
        t++;
    }
    
    //Initialize the output space (outspace)
    
    for (int i = 0; i < geosize; i++)
    {
        outspace[i][0] = geospace[i][0]; //the plant identifier
        outspace[i][1] = 0.0;
        outspace[i][2] = 0.0;
        outspace[i][3] = 0.0;
        outspace[i][4] = 0.0;
        outspace[i][5] = 0.0;
        outspace[i][6] = 0.0;
    }
    
    fclose(ifile_geodata);
    
    
    //Read in the border crossings file
    
    fscanf(ifile_lbc, "%d", &bordersize);
    double **borderspace = alloc_large_array(bordersize, 4);
    
    fprintf(stdout, "%s: Number of border crossings to be used = %d\n", argv[0], bordersize);
    
    t = 0;
    
    while (fscanf(ifile_lbc, "%lf %lf %lf %lf", &borderspace[t][0], &borderspace[t][1], &borderspace[t][2], &borderspace[t][3]) != EOF)
    {
        //We first have to convert the latitude and the longitude into radians
        borderspace[t][0] = DEG2RAD(borderspace[t][0]);
        borderspace[t][1] = DEG2RAD(borderspace[t][1]);
        
        t++;
    }
    
    fclose(ifile_lbc);
    
    
    
    //*****************************************************
    //STEP [3]: Compute the minimum distances
    //*****************************************************
    
    
    //Vector to hold the distances
    std::vector<double> temp_vect (geosize, 0.0);
    
    double mindist = 0.0;
    
    fprintf(stdout, "%s: Starting to compute the minimum distances...\n", argv[0]);
    
    FILE *ofile;
    
    if ((ofile = fopen(ofile_out, "w")) == NULL)
    {
        fprintf(stderr, "%s: ERROR - File %s could not be created, aborting.\n", argv[0], argv[2]);
        return(1);
    }
    
    for (int i = 0; i < geosize; i++)
    {
        //Reinitialize temp distance array to zero...
        for (int j = 0 ; j < geosize; j++)
        {
            temp_vect[j] = 0.0;
        }
        
        ilat = geospace[i][1];
        ilon = geospace[i][2];
        
        for (int j = 0; j < geosize; j++)
        {
            if (i != j)
            {
                jlat = geospace[j][1];
                jlon = geospace[j][2];
                
                if ((ilat != jlat) || (ilon != jlon))
                {
                    distance = (acos(sin(ilat) * sin(jlat) + cos(fabs(ilon - jlon)) * cos(ilat) * cos(jlat)) * EARTH_RADIUS);
                }
                else
                {
                    distance = 0.0;
                }
            }
            else
            {
                distance = 0.0;
            }
            
            temp_vect[j] = distance;
        }
        
        //Sort the vector in ascending order, the first x elemens contain the x smallest distances
        std::sort(temp_vect.begin(), temp_vect.end());
        
        
        mindist = 0.0;
        
        for (int j = 0; j < number1; j++)
        {
            mindist += temp_vect[j];
        }
        
        outspace[i][1] = mindist / (float)number1;
        
        
        mindist = 0.0;
        
        for (int j = 0; j < number2; j++)
        {
            mindist += temp_vect[j];
        }
        
        outspace[i][2] = mindist / (float)number2;
        
        
        mindist = 0.0;
        
        for (int j = 0; j < number3; j++)
        {
            mindist += temp_vect[j];
        }
        
        outspace[i][3] = mindist / (float)number3;
        
        
        mindist = 0.0;
        
        for (int j = 0; j < number4; j++)
        {
            mindist += temp_vect[j];
        }
        
        outspace[i][4] = mindist / (float)number4;
        
        
        mindist = 0.0;
        
        for (int j = 0; j < number5; j++)
        {
            mindist += temp_vect[j];
        }
        
        outspace[i][5] = mindist / (float)number5;
        
        
        //Min distance to US border crossing
        
        mindist = 9999999.99;
        
        for (int k = 0; k < bordersize; k++)
        {
            blat=borderspace[k][0];
            blon=borderspace[k][1];
            
            distance = acos(sin(ilat) * sin(blat) + cos(fabs(ilon - blon)) * cos(ilat) * cos(blat)) * EARTH_RADIUS;
            
            if (distance < mindist)
            {
                mindist = distance;
            }
        }
        
        outspace[i][6] = mindist;
        
        
        //Output structure: plantID [number1-number5 distances] mindistUS
        fprintf(ofile, "%lf %lf %lf %lf %lf %lf %lf\n", outspace[i][0], outspace[i][1], outspace[i][2], outspace[i][3], outspace[i][4], outspace[i][5], outspace[i][6]);
        
        if (i%100 == 0)
        {
            printf("Iteration %d done\n", i);
        }
        
    }
    
    fclose(ofile);
    
    
    //Cleanup step: Free the dynamically allocated memory
    
    geospace = free_large_array(geosize, geospace);
    outspace = free_large_array(geosize, outspace);
    borderspace = free_large_array(bordersize, borderspace);
    
    finish = clock();
    
    printf("%s: Elapsed time = %lud seconds\n", argv[0], (finish-start)/CLOCKS_PER_SEC);
    
    return(0);
}
